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CYLINDRICAL  BESSEL  FUNCTIONS  FOR  A  LARGE 
RANGE  OF  COMPLEX  ARGUMENTS 


1.  INTRODUCTION 

Many  problems  in  physics  and  engineering  require  the  use  of  so-called  special  functions;  an 
important  subset  known  as  Bessel  functions  is  essential  to  the  study  of  scattering  of  waves  by  an 
insonified  object.  This  area  of  research  is  of  particular  interest  to  those  groups  concerned  with  the  non¬ 
destructive  evaluation  of  materials,  underwater  acoustics,  and  elastic  wave  propagation. 

A  major  difficulty  arises  when  these  insonified  bodies  possess  intrinsic  absorption  properties. 
Determining  the  scattering  of  acoustic  waves  from  submerged  absorbing  bodies  requires  the  calculation 
of  complex  Bessel  functions  of  the  first  and  second  kinds,  for  orders  ranging  from  zero  to  a  value 
approximately  equal  to  the  magnitude  or  "absolute  value"  of  the  argument  [1,2].  Due  to  the  longitudi¬ 
nal  and  shear  properties  of  the  viscoelastic  materials  used,  the  magnitude  of  the  arguments  for  which 
the  Bessel  functions  are  to  be  calculated  can  be  extremely  large  [3].  To  date,  no  effective  way  has  been 
available  to  treat  Bessel  functions  of  such  large  complex  arguments  for  such  a  wide  range  of  integer 
orders. 

For  arguments  with  large  positive  or  negative  imaginary  parts,  the  values  of  the  Bessel  functions 
of  the  first  kind,  Jn(z)  ("common  Bessel  functions"),  and  second  kind,  K„(z)  ("Neumann  functions"), 
are  very  large  at  the  low  orders,  and  for  sufficiently  big  positive  or  negative  imaginary  parts  U„(z) 
approaches  Y„(z).  These  factors  create  several  problems: 

a.  when  the  values  of  the  common  Bessel  and  Neumann  functions  at  low  orders  become  larger 
and  closer  to  one  another,  matrix  singularities  arise  in  the  calculation  of  the  absorption  prop¬ 
erties,  which  lead  to  the  introduction  of  quantities  too  small  to  express  in  fixed  length  com¬ 
puter  words; 

b.  to  generate  an  accurate  table  of  Neumann  functions  by  a  recurrence  relationship,  for  succes¬ 
sive  values  of  the  order,  a  combination  of  backward  and  forward  recursion  would  have  to  be 
used;  and 

c.  the  function  values  can  become  so  large  as  to  exceed  the  permissible  floating-point  exponent 
range  of  the  computer. 

The  next  section  discusses  the  steps  the  author  has  taken  to  overcome  these  problems. 

2.  APPROACH 

To  circumvent  the  first  two  problems,  Bessel  functions  of  the  third  kind,  Hj'Hz),  /  —  1,  2 
("Hankel  functions"),  rather  than  Neumann  functions,  are  calculated  when  ris  large.  The  values  of  the 
Hankel  functions  are  small  in  magnitude  at  low  orders,  thus  reducing  the  chance  of  a  matrix  singularity; 
moreover,  for  orders  of  value  greater  than  1  they  can  be  generated  by  the  forward  relationship  alone. 

As  implemented,  the  quadrature  method  (Simpson's  Rule)  used  to  calculate  the  Hankel  functions 
does  not  yield  accurate  results  for  small  arguments;  therefore  the  absorption  program  is  written  to  use 
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Neumann  functions  for  small  arguments  and  Hankel  functions  for  large  arguments.  A  general-purpose 
subroutine  was  written,  however,  that  can  generate  all  three  functions  for  any  complex  argument  and 
any  positive  integer  order.  If  Y„(z)  for  0  <  n  ^  Wand  large  |z|  is  to  be  generated,  the  subroutine  first 
calculates  J„(z)  and  either  //„(l)(z)  or  //„<2)(z).  Then  it  uses  the  relationship: 

tf„(l)(z)  -  J„(z)  +  iYJz ) 


or 

//„<2)(z)  -  Jn(z )  -  /T„(z) 

to  generate  Y„(z).  (In  the  program  H^Hz)  is  calculated  when  z  is  in  the  first  or  second  quadrant; 
H„a)(z)  is  calculated  when  z  is  in  the  third  or  fourth  quadrant.)  Alternately,  if  |z  I  is  small  and  f/„(,)(z) 
is  to  be  generated,  Jn(z )  and  Y„(z )  are  calculated  first  and  then  one  of  these  two  relationships  is  used 
to  solve  for  H^Hz). 

The  subroutine  is  designed  so  that  a  table  of  values  for  each  of  the  functions  is  generated,  for 
successive  integer  values  of  n  from  zero  to  the  order  specified.  Generally,  for  functions  of  order 
n  >  1,  answers  are  generated  by  using  the  recursion  relationship  [4,5]: 

C„_i(z)  +  Cv+X(z)  -  (2 viz)  C„ (z) . 

The  one  exception  occurs  for  Jj(z)  to  Jn(z )  when  |z|  <  0.5  and  the  maximum  order  <5;  see  Part 
3b(l),  Methods  of  Computation.  For  small  arguments,  entire  tables  of  both  Jn(z)  and  Y„(z)  are  calcu¬ 
lated;  then  values  of  //„(/)(z)  are  generated  from  these.  For  large  arguments,  tables  of  Jn(z )  and 
Hj'Hz)  are  calculated;  and  then  the  Y„(z)  table  is  generated  by  using  them.  In  this  way,  one  avoids 
the  not  so  simple  problem  of  generating  an  accurate  table  of  Y„(z)  for  extremely  large  arguments  by 
the  recursion  operation.  The  difficulty  encountered  if  one  were  to  attempt  to  use  the  recursion  rela¬ 
tionship  to  compute  a  Yn(z)  table  arises  because  the  recursion  operation  works  best  if  each  successive 
function  value  is  greater  in  magnitude  than  the  one  preceding.  And  |  Y„  (ziarge)  I  is  initially  a  decreasing 
function,  but  at  a  value  of  n  roughly  equal  to  {z|,  it  begins  to  increase.  (A  similar  difficulty  arises  in 
the  recurrence  calculation  for  the  spherical  Bessel  function,  y„(z)  [6].)  This  means  that  the  Y„(z)  table 
would  have  to  be  calculated  backward,  roughly,  from  n  —  |z|  to  n  —  0,  and  forward  from  n  —  \z\  to 
maximum  /r,  and  therefore  one  would  have  to  be  able  to  predict,  fairly  accurately,  where  the  minimum 
absolute  value  of  the  function  occurs.  On  the  other  hand,  |/„(z)|,  I  TB(zsnM„)|,  and  |//„(,)(z)|  are  either 
monotonically  increasing  or  monotonically  decreasing  functions,  so  the  recursion  operation  can  begin  at 
n  —  1  (or  n  —  n^).  See  Figs.  1  and  2. 

Figures  1(a),  (b),  and  (c)  show  the  normalized  curves  of  the  real  part  for  each  of  the  three  kinds 
of  Bessel  functions,  for  an  argument  equal  to  201  + 150/  and  a  range  of  integer  orders  from  zero  to  462. 
Figure  1(d)  shows  the  three  curves  superimposed.  Unfortunately,  if  the  maximum  order,  n,  were 
made  larger  than  462  to  show  better  the  oscillations  in  Y„(z)  and  tf„(1)(z)  at  large  orders,  then  the 
Y„(z)  values  at  low  orders  would  appear  to  vanish  because  the  function  values  at  n  >  462  are  so  much 
larger,  relatively.  Normalized  curves  of  the  imaginary  part  follow  the  same  general  pattern. 

Figure  2  gives  another  illustration  of  the  behavior  of  these  functions.  Here  log  |/„(z)|, 
log  |K„(z)l,  and  log  |//„(1)(z)|  are  plotted  for  six  different  values  of  the  argument,  z.  Note  that  the 
turning  points  for  the  Neumann  functions,  in  Fig.  2(b),  occur  at  values  close  to  but  always  smaller  than 
zero,  i.e.  they  have  a  function  value  less  than  one.  Also  note  that  when  the  three  sets  of  plots  are 
superimposed,  as  in  Fig.  2(d),  in  this  scale  the  Hankel  curves  appear  to  be  a  reflection  of  the  common 
Bessel  curves,  and  the  Neumann  curves  coincide  with  the  common  Bessel  curves  before  the  Neumann 
turning  point  and  with  the  Hankel  curves  after  the  turning  point. 

The  numerical  problem  of  exceeding  the  exponent  range  of  the  computer  has  been  solved  by  a 
programmed  scaling  of  each  step  of  the  calculation,  where  necessary,  over  and  above  the  built-in  digital 
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Fig.  1  -  Computer-generated  plots  (normalized),  Z  -  201  +  150/.  (a)  Bessel  functions;  (b)  Neumann  functions;, 
(c)  Hankel  functions;  (d)  Bessel,  Neumann,  and  Hankel  functions  superimposed. 
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Fig.  2  —  Compuier-generaied  plots  (logarithmic),  (a)  Bessel  functions;  (b)  Neumann  functions;  (c)  Hankel  functions; 

(d)  Bessel,  Neumann,  and  Hankel  functions  superimposed. 

scaling  that  is  effected  automatically  by  the  machine's  double-precision  floating-point  arithmetic.  In  our 
case,  the  values  of  the  functions  are  scaled  by  multiples  of  10±7°,  but  the  size  of  this  factor  depends 
really  upon  the  exponent  range  of  the  particular  computer  in  use. 

Theoretically  there  should  be  no  limits  on  the  values  of  n,  x,  or  y,  but  at  the  present  time  the 
subroutines  have  been  tested  only  for  integer  n  on  the  interval  0  <  n  <  3010,  and  on  the  rectangular 
region  defined  by  the  complex  points  ±3000±  3000/. 

3.  METHODS  OF  COMPUTATION 

a.  Procedures 

The  necessity  for  using  various  methods  or  a  combination  of  methods  was  indicated,  as  there 
appeared  to  be  no  single  superior  one  [7,8].  Figure  3  shows  the  number  of  different  procedures  used  to 
generate  the  three  function  tables  and  the  portion  of  the  complex  plane  in  which  each  of  these  methods 
is  applied.  The  precise  outline  of  these  procedures,  including  the  formulas  used,  is  given  in  Part  3b. 
Notice  that  when  5.0  <  [y|  <  10.0,  there  is  a  choice  of  one  of  three  possible  procedures.  The  integral 
method,  given  in  Part  3b(6),  is  always  employed  if  the  Hankel  function  is  to  be  calculated;  otherwise, 
the  power  series,  Part  3b(2),  or  the  asymptotic  expansion,  Part  3b(3),  is  used. 
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Fig-  3  —  Procedural  regions  in  the  Z-plane  (see  Part  3b  for  definitions  of 
(I),  (2),  (3),  and  (6)) 


b.  Definitions 

(1)  If  \z |  <  0.5  and  maximum  order  <5;  /0(z),  J\{z),  and  Y\(z)  are  calculated  by 

means  of  a  power  series  [9,10]: 

u0(p,<f>)  -  £  (-D*  y  cos  2k * 

(P.4>)  -  I  (~D*  sin  2k<f> 

d 

uy(p,(f>)  -  -  —  [«0  cos  tf>  +  v0  sin  0] 
op 

i>i(p,0)  —  —  [«o  sin  <f>  -  v0  cos  ^] 

Op 

Jr.ipe1*)  -  w„(p,<£)  +  iv„  (p,rf>) 


U0(p,<f>)  -  |mo(p.0)|>'  +  In  -  <^*/0(p.<A)|  +  5o(p,0) 

K0(p,^)  -  L(p,*)  fy  +  In  ^-1  +  +  T0(p,<t> ) 
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l«l(p,0) 

|y  +  ln^- 

|  -  0^l(p,</l)j 

y  +  In  j’j 

+  4>U\(p,<f>) 

So<p.*)  -  ±  (l+-j  +  ...+  -j  cos  2  ** 

r°<p-*) '  »  Ji  (1  +  T++i) si"  2‘* 


2|l  +  y  +  .. 

+  1 

+  -i— 

k 

k+l 

2  1  +  1  +  •• 

,  +  I 

+  \ 

1  2 

k 

k+l 

Y„(pe14)  -  £/„(p,$)  + 

4(2)  to  J„(z)  are  calculated  from:’ 

y  (2)  „  i.  y  _ (— z2/4)* 

"  2]  is  *!T(»  +  it  +  l)- 

F2(2)  to  Yn{z)  are  calculated  by  using  the  forward  recursion  relationship: 

Wz)-  (2»/z)f,,(z)-  y„_,(z). 

For  Hq*  (z)  to  //„(,)(z)  one  of  the  following  relationships  is  used: 

//„(l)(z)  -  y„(z)  +  /r„(z) 
or 

h„(2Hz)  -  y„(z)  -  /r„(z). 

(2)  If  Izf  ^  0.5,  |x I  <  16.0,  and  |.y|  <  5.0  or  if  |z I  <  0.5  and  maximum  order  >5;  J0(z), 
Jyiz),  Y0(z),  and  f((z)  are  calculated  by  means  of  the  power  series  in  Part  3b(l). 

J2(z)  to  J„(z)  are  calculated  by  using  the  backward  recursion  relationship  [11-13]: 

i(z)  -  (2 n/z)J„(z)  -  /„+|(z). 

Y2(z)  to  Y„(z)  are  calculated  by  forward  recursion,  as  in  Part  3b(l).  //0(,)  (z)  to  H}n(z) 
are  calculated  by  the  same  method  as  in  Part  3b(l). 

(3)  If  |x|  >  16.0  and  [yl  <  5.0;  /0(z),  /j(z),  fo^),  and  Y\(z)  are  calculated  by  means  of 
an  asymptotic  expansion  [9,10]: 

1  .  1 1/2  it  1  f  it 


4(2)  ~  |f„(z)  cos  2  -  ~  -  ~  -  Qv(z)  sin  2  -  -  —■  | 

Yv(z)  —  |-^J  1^(2)  sin  z  -  ~  -  y  +  Q„(z)  cos  z  -  ~  ~  J 


la 

•rl m  --  «  #  •  *  ‘  m  ■  m  m  •  ■  m  • 

1  ;  • » •< 

.  •  a  W  g  *  *  •  #  *  >  ’  »  «  «  « 

■  <Y  -V  .A 
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where 


Pjz)  - 

CM*)  - 


1  +  I 

A-l 


(—  \)k(4v2  —  12)(4i/2  —  32)  ...  (4v2  —  {4A:  —  l}2) 


I 

*- 1 


(2fc)!26*z2* 

(— 1)*-m  (4^2  _  | Z) (4^2  —  3^)  (4^2  _  {4^  —  3)7) 

(2k  —  1)!26*-3z2*-1 


J2(z)  to  J„(z)  are  found  by  backward  recursion;  Y2(z)  to  Y„(z)  are  found  by  forward 
recursion.  ( z )  to  A/„(f)(z)  are  calculated  by  using  the  same  method  as  in  Part  3b(l). 

(4)  If  |x|  <  16.0,  5.0  <  |y|  <  10.0,  and  the  Hankel  functions  are  not  needed,  the  common 
Bessel  and  Neumann  functions  are  calculated  as  in  Part  3b(2).  If,  on  the  other  hand,  the 
Hankel  function  is  to  be  computed,  the  common  Bessel,  Neumann,  and  Hankel  functions 
are  calculated  by  using  the  relationships  given  in  Part  3b(6),  q.v. 

(5)  If  |x|  >  16.0,  5.0  <  |j/|  <  10.0,  and  the  Hankel  functions  are  not  needed,  the  common 
Bessel  and  Neumann  functions  are  calculated  as  in  Part  3b(3).  If,  however,  the  Hankel 
function  is  to  be  computed,  the  common  Bessel,  Neumann,  and  Hankel  functions  are  cal¬ 
culated  by  using  the  relationships  given  in  Part  3b(6),  q.v. 

(6)  If  M  >  10.0  (or  greater  than  5.0,  if  Hankel  functions  are  needed),  J0(z)  and  J\(z)  are 
calculated  by: 

1  cv 

J„(z)  -  —  I  cos  ( z  sin  9  -  nO)  d9. 

it  Jo 


For  (z)  and  ff{°  (z),  where  y  is  positive: 


H„0)(z) 


-2  ie~'”v  (z/2)”  r” 

Vw  r (1/2  +  n) 


,2i.e«r(r2+l)'/J 

(t2  +1)'/2 


and  where  y  is  negative: 


H™(z) 


2/(-l)"  (z/2)”  r” 

>/ir  r (1/2  +  n)  Jo 


/2de-fe(/2+  1)1/J 

(t 2  + 1),/2 


For  each  of  the  integrals,  Simpson's  quadrature  method  is  used;  the  number  of  subdivi¬ 
sions  necessary  is  determined  by  convergence  criteria. 

J2(z)  to  J„(z)  are  found  by  backward  recursion;  Hjn  (z)  to  Hf‘Hz)  are  found  by  forward 
recursion.  For  Y0(z)  to  Y„(z),  one  of  the  following  relationships  is  used: 

Yn(z)  -  Un(z)  -  i//„(1)(z) 

or 

Y„(z)  -  -  U„(z)  +  iH™(z). 


4.  VERIFICATION 

Because  published  tables  often  a?  limited  ir  ige  or  lack  precision,  the  accuracy  of  the  derived 
function  values  was  checked  with  the  ian  ationship  (4]: 
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Wz)  y„(z)  -  J„(z)  Yn+i(z)  =  2/ (irz) 

[yn(z)  //ji*!  (z)  -  y„  +  1(z)  //„(1)(z)]/  =  2/(?rz) 

yn+1(z)  r„(z)  -  yn(z)  k„+1(z)  -  2/(ttz) 

-u„(z)  //„(2)i  (z)  -  J„+i(z)  //„(2>(z)l/  -  2/ (nz). 


Table  1  gives  an  example  of  the  calculated  values  of  y„(z),  T„(z),  and  //„(2)(z),  to  16  digits,  and  of  the 
corresponding  Wronskians  for  an  argument  z  -  3000-3000;  and  orders  n  =■=  0,  1,  2,  ,  9.  The 

number  in  the  last  column  is  the  value  of  u  in  the  expression  1070“  that  is  used  to  multiply  each  value 
on  that  line.  For  example,  the  fourth  line  states  that 


Table  1  —  Example  of  the  Calculated  Values  of  Jn(z),  Y„(z),  and  /f„(2)(z),  to 
16  Digits,  and  of  the  Corresponding  Wronskians  for  an  Argument  z  =  3000  - 
3000i  and  Orders  n  —  0,  1,  2 . 9 

l  -  0.30000000000000000  04  -0.30000000000000000  04 

2/CPI *Z)  =  0. 1061032  9539459690-03  0  0.1 061 0 3295 3945 9690-0 3  0 

REAL  PART  IMAGINARY  PART 


N 

= 

0  JCZ) 

= 

-0.38286469325356710  41 

18 

0.26970776985386840  41 

18 

Y<2> 

= 

0.26970776985386840  41 

18 

0. 38286469325  35 67 ID  41 

18 

H(Z) 

= 

-0.27370745790819070-45 

-18 

-0.15784623445389830-44 

-18 

WRONSKIAN 

= 

0.10610329539459630-03 

0 

0.1061 032 953947 169D-0 3 

0 

N 

= 

1  JCZ) 

= 

0.26965338615636280  41 

18 

0. 38 2855265 38 36 766D  41 

18 

YCZ) 

S 

0.382855265 38 36766 D  41 

18 

-0.26965338615636280  41 

IE 

HC2  ) 

s 

0.15785710760114310-44 

-18 

-0.27 3861 794431 5793 D-45 

-IE 

WRONSKIAN 

3 

0. 1 061 0  3295  394  596  30-0  3 

0 

0.10610329539471690-03 

C 

N 

3 

2  JCZ) 

a 

0.38282695929382470  41 

18 

— 0.26 94 90266 97002 19D  41 

IE 

YCZ) 

3 

-0.26949026697002190  41 

18 

-0.38282695929382470  41 

IE 

HCZ) 

3 

0.27432493553167170-45 

-18 

0.15788972476328430-44 

-IS 

WRONSKIAN 

3 

0.10610329539459630-03 

0 

0.1061Q32953947169D-03 

C 

N 

= 

3  JCZ) 

* 

-0. 2692 1850 80 055 2 02 D  41 

18 

-0.382 779707588794 ID  41 

IE 

YCZ) 

= 

-0.3827797075887941D  41 

18 

0.26 921 85 08005 5 202D  41 

18 

HCZ) 

3 

-0.15794407908861650-44 

-18 

0.2750972758870223D-45 

-IE 

WRONSKIAN 

= 

0.10610329539459630-03 

0 

0.10610329539471690-03 

C 

N 

= 

4  JCZ) 

3 

-0.38271339809424140  41 

18 

0.26883826875442760  41 

IE 

YCZ) 

3 

0.26883826875442760  41 

18 

0.38271339809424140  41 

IE 

HCZ) 

X 

-0.2761 7947359844490-45 

-18 

-0.15802015911478420-44 

-IE 

WRONSKIAN 

3 

0.10610329539459630-03 

0 

0. 1061032 953947169D-03 

C 

N 

3 

5  JCZ) 

* 

0.26834977244972210  41 

18 

0.38262787408300770  41 

It 

YCZ) 

3 

0.38262787408300770  41 

18 

-0.26834977244972210  41 

18 

HCZ) 

3 

0.15811794870428970-44 

-18 

-0.27757245064001730-45 

-IE 

WRONSKIAN 

3 

0.10610329539459630-03 

0 

0.10610329539471690-03 

C 

N 

3 

6  JCZ) 

3 

0.38252293459151940  41 

18 

-0.26775330601020640  41 

IE 

YCZ) 

3 

-0.26775330601020640  41 

18 

-0.38252293459151940  41 

18 

HCZ) 

3 

0.27927739349458310-45 

-18 

0.15823742695418470-44 

-IS 

WRONSKIAN 

3 

0.10610329539459630-03 

0 

0.10610329539471690-03 

0 

N 

3 

7  JCZ) 

3 

-0.26704921996851870  41 

18 

-0.38239833482584510  41 

18 

YCZ) 

* 

-0.382398 3348 25845 ID  41 

18 

0.26704921996851870  41 

18 

HCZ) 

3 

-0.15837856807949920-44 

-18 

0.28129575396609020-45 

-18 

WRONSKIAN 

3 

0.10610329539459630-03 

0 

0.10610329539471690-03 

C 

N 

* 

8  JCZ) 

X 

-0.38225378665685230  41 

18 

0.26623792838235290  41 

18 

YCZ) 

X 

0.26623792838235290  41 

18 

0.3822537866568523D  41 

18 

HCZ) 

3 

-0.28362925017569230-45 

-18 

-0.15854134127044480-44 

-18 

WRONSKIAN 

a 

0.10610329539459630-03 

0 

0.10610329539471690-03 

C 

N 

3 

9  JCZ) 

3 

0.2653 1990 8728414 iO  41 

18 

0.38208895920377980  41 

18 

YCZ) 

3 

0.38208895920377980  41 

18 

-0.26531990872841410  41 

18 

HCZ) 

3 

0* 15872 57 10522 8402 D- 44 

-18 

-0.28627986773377050-45 

-18 

WRONSKIAN 

3 

0.10610329539459630-03 

0 

0.10610329539471690-03 

C 
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y0(3000-3000/)  =  (-0.382864  ...  x  1041  +  0.269707  ...  x  10 4,i)  x  10,70xl8\ 
in  other  words, 

J0(3000-3000i)  =  -0.382864  . . .  x  101301  +  0.269707  . . .  x  101301/. 

The  seventh  line  gives  the  Wronskian  for  the  calculated  function  values.  This  should  be  compared  with 
the  value  of  2/(ttz)  given  in  the  second  line.  Table  2  lists  the  function  values  and  Wronskians  for  the 
same  argument,  but  for  orders  n  =  3000,  3001,  3002 .  3009. 

Although  it  would  be  impractical  to  describe  here  all  the  tests  that  were  performed  to  determine 
the  best  range  for  each  of  the  procedures  covered  in  Part  3,  some  discussion  might  be  of  interest.  One 
set  of  tests  consisted  of  generating  J„(z)  and  Y„(z)  (or  J„(z),  Tn(r),  and  H„(z))  for 
0  <  n  <  24  near  the  points  of  procedure  change,  e.g.  z  =  0  ±  5i,  |z  |  =  0.5,  z  =  ±16  ±  5i,  etc. 
(see  Fig.  3).  The  computed  function  values  and  their  corresponding  Wronskians  were  examined  at  n  = 
24.  Whether  comparing  the  function  values  generated  by  the  integral  method  to  those  generated  by  the 
series  method  or  to  those  generated  by  the  asymptotic  method,  the  .^(zFs  always  agreed  to  at  least  12 
figures  and  frequently  to  16  figures.  And  for  Y2^z)  agreement  ranged  from  4  figures  at  z  =  15.95  + 
lO.Oi  to  12  figures  at  z  =  0.01  +  5.01  i.  The  accuracy  of  the  Wronskians  varied  from  a  low  of  4  figures 
for  z  =  15.95  +  lO.Oi,  if  the  Bessel  and  Neumann  functions  were  used,  to  a  high  of  15  figures  for 
several  values  of  z.  Interestingly,  when  the  routine  was  temporarily  modified  in  such  a  way  that 
K0(l 5.95  +  lO.Oi)  was  calculated  by  the  series,  asymptotic,  and  integral  methods,  all  three  calculations 
agreed  to  14  figures  for  n  =  0,  and  P0(z)P\(z)  +  Qo(z)Qx(.z)  equaled  1.0  to  14  figures.  The  three 
methods  gave  values  that  agreed  to  7  figures  for  n  =  20.  This  indicates  that  the  calculation  of  T0(15.95 
+  lO.Oi)  by  the  series  method  is  sufficiently  accurate  but  that  as  n  increases  the  forward  recursion  loses 
accuracy. 

The  P„(z)P„+\(.z)  +  Qn(z)Qn+l(z )  —  1.0  test  was  also  made  using  the  P  and  Q  results  for 
z  =  150.0  +  4.95i  at  n  —  0. 

The  test  answer  equals  0.1000000000000000  x  10'  -  0.6286572655403010  x  10-22i. 

5.  PORTABILITY 

The  results  shown  in  Tables  1  and  2  were  generated  by  a  subroutine  written  in  double-precision 
FORTRAN  and  run  on  the  Texas  Instruments,  Inc.  computer,  ASC.  This  computer  has  a  word  length 
of  32  bits,  with  an  exponent  range  of  approximately  ±75.  The  subroutine  was  easily  modified  for  run¬ 
ning  on  a  Digital  Equipment  Corp.  PDP-11.  This  computer  also  has  a  word  length  of  32  bits;  however, 
because  the  exponent  range  is  approximately  ±36,  the  scale  factor  had  to  be  reduced.  Even  though  the 
PDP-11  does  not  provide  double-precision  complex  arithmetic,  no  accuracy  was  lost  since  only  real 
arithmetic  is  employed  in  the  subroutine.  There  is  less  internal  storage  space  in  the  PDP-11,  so  the 
maximum  allowable  order  and  argument  sizes  are  approximately  r.  <  230  and  |z|  <  230. 

The  integral  method  is  always  employed  if  the  Hankel  function  is  to  be  calculated  prior  to  the 
Neumann  function.  This  approach  was  taken  because,  although  the  asymptotic  expansion  portion  can 
be  modified  rather  easily  to  give  H0(z)  and  H\(z)  (in  fact,  instructions  for  doing  so  are  incorporated  in 
the  subroutine),  scaling  here  would  have  presented  many  more  problems  than  scaling  the  integral 
method.  If  one  has  access  to  a  computer  with  a  very  large  exponent  range,  e.g.  the  DEC  VAX-750, 
the  asymptotic  expansion  could  be  used  for  a  much  greater  range  of  z-values  than  shown  in  Fig.  3. 

6.  FUTURE  USE 

Up  to  the  present  time,  these  Bessel  function  subroutines  have  been  used  in  the  determination  of 
the  acoustic  reflection  from  absorbing  infinite  cylinders.  These  subroutines  will  become  even  more 
valuable  as  the  work  being  done  on  finite  cylinders  is  expanded  to  include  the  absorption  properties. 
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Table  2  —  Example  of  the  Calculated  Values  of  J„(z),  Y„(.z),  and  HnQ}(z).  to 
16  Digits,  and  of  the  Corresponding  Wronskians  for  an  Argument  z  =  3000  - 
3000i  and  Orders  n  =  3000,  3001,  3002 .  3009 

2  =  0.30000000000000000  04  -0.30000000000000000  04 

2/CPl*Z)  -  0. 106103295 39459690—03  0  0.10610329539459690-03  0 


N  =  3000  JCZ)  = 

T(Z)  = 
HCZ)  = 
URONSKIAN  = 
N  =  3001  JCZ)  = 

TCZ)  = 
HCZ)  = 
URONSKIAN  = 
N  =  3002  JCZ)  * 

rcz>  = 
HCZ)  = 
URONSKIAN  * 
N  =  3003  JCZ)  = 

TCZ)  = 
HCZ)  * 
URONSKIAN  - 
N  *  3004  JCZ)  = 

TCZ)  * 
HCZ)  = 
URONSKIAN  = 
N  *  3005  JCZ)  = 

TCZ)  = 
HCZ)  « 
URONSKIAN  = 
N  =  3006  JCZ)  = 

TCZ)  = 
HCZ)  = 
URONSKIAN  * 
N  »  3007  JCZ)  * 

TCZ)  * 
HCZ)  = 
URONSKIAN  a 
N  *  3008  JCZ)  = 

TCZ)  * 
HCZ)  * 
URONSKIAN  a 
N  a  3009  JCZ)  a 
TCZ)  * 
HCZ)  a 
URONSKIAN  * 
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REAL  PARI 

-0.40619757734038320  53  13 

0.51063286538286730  54  13 

0.11168513570573600-57  -13 
0.10610329539459700-03  0 

0.25967407715839060  54  13 

0.1527705590799730D  54  13 

0.20 82 7660 7 985 67 45 0-57  -13 
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0. 1475589103184821 D  54  13 

-0.98050747599090880  53  13 

-0.13285865991062630-58  -13 
0.10610329539459700-03  0 

-0.13900679468872620  53  13 

-0.10322939091876890  54  13 
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0.10610329539459710-03  0 

-0.58140870157 135 86 D  53  13 

-0.19196452858938300  53  13 

-0.88697492758311060-57  -13 
0.10610329539459710-03  0 

-0.25095 663 7 19055 85 D  53  13 

0.25788951472006680  53  13 

0.48025666110424300-57  -13 
0.10610329539459710-03  0 

0.71714472740882430  52  13 

0.19890896091477370  53  13 

0.32830139436481110-56  -13 
0.10610329539459700-03  0 

0.1235077600403195D  53  13 

0.13275165802900650  52  13 

0.34992736822757800-56  -13 
0.10610329539459700-03  0 

0.38775330883090380  52  13 

-0.61806874911252740  52  13 

-0.42975271491301330-56  -13 
0.10610329539459700-03  0 

-0.22 65 73350 305 24 SOD  52  13 
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0.40619757734038320  53  13 

-0.81937857861700540-58  -13 
0.10610329539471770-03  0 

0.15277055907997300  54  13 

-0.25967407715839060  54  13 

0. 10991012T09796370-5T  -13 
0.10610329539471770-03  0 

-0.98050747599090880  53  13 

-0.14755891031848210  54  13 
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0.10610329539471770-03  0 

-0.10322939091876890  54  13 

0.13900679468872620  53  13 
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0.10610329539471770-03  0 

-0.19196452858938300  S3  13 
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-0.74535161292567210-57  -13 
0.10610329539471770-03  0 

0.25788951472006680  53  13 

0.25095663719055850  S3  13 
— 0. 191 1795601 190267D-56  -13 
0.10610329539471770-03  0 

0.19890896091477370  53  13 

-0.71714472740882430  52  13 

-0. 68857322539382860-57  -13 
0.10610329539471770-03  0 

0.13275165802900650  52  13 

-0.1 23507 760040 31 95D  53  13 

0.45114252008810570-56  -13 
0.10610329539471770-03  0 

-0.61806874911252740  52  13 

-0.38775330883090380  52  13 

0.87179637392780290-56  -13 
0.10610329539471770-03  0 

-0.36368127281804770  52  13 

0.22657335030524800  52  13 

-0.79200779826099200-58  -13 
0.10610329539471770-03  0 
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